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The Karlqvist approximation signaling the historical beginning of magnetic recording head the¬ 
ory is reviewed and compared to various approaches progressing from Green, Fourier, Conformal 
mapping that obeys the Sommerfeld edge condition at angular points and leads to exact results. 
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I. INTRODUCTION 

Karlqvist seminal historical paper® is an important 
landmark and the first building block of magnetic record¬ 
ing theory. Despite all ensuing developments in recording 
theory it is still used as a standard reference with respect 
to any model aiming at the description of recording or 
replay head fields. 

Karlqvist work in based on magnetostatics and Green’s 
function theory of partial differential equations since the 
basic description of head field physics is based on Laplace 
or Poisson equations the magnetostatic potential must 
satisfy under certain conditions. 

Magnetic playback theory deals with time dependence 
since tapes and hard disk platters move with respect to 
head in contrast with our work that deals with static 
head description. 

We concentrate on the description of magnetostatic in¬ 
teractions in 3D simplifying them to 2D in order to fully 
understand the Karlqvist approximation and its many 
aspects that lie at the roots of recording theory. 

This work can be taught as an application chapter 
in a general Electromagnetism course at the Graduate 
level or in an undergraduate class of Partial Differential 
Equations since physicists might be interested in apply¬ 
ing some mathematical methods to solve Laplace and 
Poisson equations or understanding issues surrounding 
the development of magnetic recording advances in hard 
disks and storage media. 

Karlqvist approximation is based on a simplified de¬ 
scription of the recording head physics (see fig. and 
fig-d to 2D such that it consists of two semi-infinite 
poles made from a very soft material possessing an in¬ 
finite permeability. The poles are separated by a gap 
region and face a semi-infinite plane devoid of any mag¬ 
netic charges. 

In addition a number of simplifying assumptions are 
made in order to obtain a fully analytic solution to the 
field created by the poles, such as: 

1. Uniform magnetization above poles leading to a 
zero volume source density. 

2. Spatially linear surface source density in the gap. 

3. Uniform field Hg in the gap all the way up to y = 0. 


In reality, it can be uniform only in the case y << 0 
far from the dFl boundary. 

4. Ignorance of Sommerfeld radiation conditiorP near 
angular boundaries such as regions around points 
X = ±a, y = 0. 

In this work we review briefly how magnetostatic the¬ 
ory is used for the description of fields created by the 
head and its interaction with recording media. It is or¬ 
ganised as follows. Section 2 describes Green’s function 
theory in 3D and its application to 2D with the Karlqvist 
approximation delivering the potential and field created 
by the recording head in the plane on the basis of a sim¬ 
plified Boundary Value Problem (BVP) described above. 
In section 3, Karlqvist solution is revised and compared 
to exact Fourier methods meant to solve the exact BVP. 
In section 4 we use conformal mapping methods to re¬ 
trieve the solution for a more realistic head accounting 
for edge condition at angular points and leading to exact 
results. Finally section 5 bears conclusions and perspec¬ 
tives of the work. 


II. GREEN’S FUNCTION APPROACH 

In the 3D magnetostatic approximation. Maxwell 
equations for fields E and H reduce to: 

Vx£; = 0, Vxff = 0 (1) 

In the absence of magnetization and poles (V -3 = 0) 
we infer existence of a magnetic scalar potential (/'(r) at 
any point r leading to H definition and Laplace equation: 

H = -V4>{r), A<()(r)=0 (2) 

In presence of a magnetization field M{r), the low¬ 
est approximation for the ’’far-held” potential (j>{r) in a 
multipole expansion is dipolar leading to a volume source 
V • M and a surface density source n ■ M where n is the 
outward normal to the hnite surface sample embodying 
M{r). 

The potential can be obtained, from volume 

source density p{r) = —V • M{r) and surface source 
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[ [(l){r)AG-GA(l){r)]dV = i [(l){r)VG-GV(j){r)]-ndS 
Jv JdV 

( 5 ) 

derived from coupling Laplace equation and the defi¬ 
nition of the Green’s function solution of: 


AG{r, r') = (5(r — r') (6) 

where S(r) is Dirac delta function. 

The radial Green’s function expression is obtained by 
solving the 3D Laplacian in spherical coordinates: 


l±f,2dG 

dr \ dr 


= 0 for r > 0 


( 7 ) 


Thus the full Green function G{r,r') = 
obtained (the angular factor 47r originates from solid an¬ 
gle subtending all 3D space) and exploited in expressing 
the potential as: 


(j)(r) = f G{r,r')p{r')dV'+(1) VG{r,r')-n'(j)s{r')dS' 

Jv JdV 

( 8 ) 

Using eq. E eq. and the condition V())(r) = 0 on dV 
we retrieve the general solution for the potential eq. 


A. Karlqvist 2D magnetostatic solution 


FIG. 1: (a) General structure of a magnetic recording head. 
Pole and core are generally made from very soft material pos¬ 
sessing a very large permeability fj,. They are excited by a 
coil that will create a magnetic field in the gap and around 
the pole surfaces, (b) 3D perspective view of a simplified 
magnetic recording head tip. 


density cr(r) = n ■ M{r) as (considering free-space per¬ 
meability fiQ = 1): 




J_ r p{r')dV' J_ r a{r')dS' 
dTT Jy \r - r'\ An Jgy \r - r'\ 


( 3 ) 


Karlqvist approximation is based on approximating 
the magnetic recording head as consisting of the tips only 
(see fig{T) with the assumption the poles are large and 
thick and made with an infinitely soft material (perme¬ 
ability /r « oo). 

A flat 2D model for the recording head displayed in 
%[2] can be made since poles are considered as infinitely 
large (along x direction as in fig. and thick (along z 
direction, see fig. separated by the gap region. 

With this assumption, one cannot use the potential 
result previously given in eq. since this is valid strictly 
in the 3D case. 

Thus we move to a fully 2D approach to solve the BVP 
posed by the magnetic head stated mathematically as: 


This result can in fact be reinterpreted as a BVP with a 
Poisson equation subjected to presence of volume source 
and surface sources: 

A(j){r) = V • AT, in V, (p(r) = n ■ M on dV (4) 

It is interesting to note that the surface density is 
equivalent to a surface potential = n ■ M{r). 

The BVP eq. can be solved with Green second iden¬ 
tity: 


A(j){r) = 0 in D, 4>ir) = <(>s(t’) on dfl (BVPl) 

( 9 ) 

This means the magnetization held is uniform yielding 
zero volume source, however we have a non-zero surface 
potential 0s(r) on 9D. 

The Green’s function of the 2D Laplace equation is 
G(r, r') = ^ln|r — r'\ obtained by solving the 2D 
Laplace equation in polar coordinates for the radial part 
(as in the 3D case, 27r corresponds to angle subtending 
the entire plane): 
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FIG. 2: (Color on-line) 2D approximation of a magnetic 
recording head. BVPl eq. |^is defined by domain fl (y > 0) 
bounded by dQ (in dashed red), the entire x axis. The poles 
are made from very soft material, H = 0 for any point on 
dQ and inside poles since Hg is finite in the gap region. Note 
vertices x = ±a, y = 0 that should be subject to Sommerfeld 
radiation condition. 


As in the 3D case, the solution to the BVP is obtained 
from Green second identity over the domains Ql and dOl 
corresponding to the previous 3D volume and surface V 
and dV respectively using an adequate 2D Green func¬ 
tion. 

The 2D Green function G{r, r') = ^ In |r — r'| is for 
the entire plane whereas in our case, the domain Ql is 
limited to the upper half-plane, thus the Green function 
should be G(r, r') = ^ In |r — r'| — ^ In \r — r^\ using 
the method of imaged, r' = {x',—y') is the image of 
r' = {x',y') in the lower half-plane. 

The explicit form of the Green function in the upper 
half-plane becomes: 


G{r, r') 



^J{x-x'f + {y-y'f 

^J{x-x’f + {y + y'f 


{x-x'f + {y-y'f 

47r (a; — x’)^ + {y + y')^ 


( 11 ) 


Using eq. 1^ eq. |^and V4>{r) = 0 on dQ: 


(j){r)= [ (l)s{r')V'G{r,r')-n'dS' (12) 

Jan 

The quantity V'G • n' over the boundary dQ is ob¬ 
tained with the outward normal n' = (0, —1) as: 


V'G • n' = 


yielding: 


dG{r, r') 


dy' 


v'=o 


y 1 

7’" [(x — x')^ + y^] 


(13) 


(j){r) = (t){x,y) 


y 

TT 



4)s{x')dx' 

{x — x')^ + y^ 


(14) 


The above result called the Poisson integral for the 
half-plane may be derived from the fact an analytic func¬ 
tion can be obtained at any point from the values it takes 
over a boundary (see for instance Ablowitz et al. 1^. It 
can be viewed as a convolution between the source func¬ 
tion (j)s{r) and the Green function. The result is retrieved 
with Fourier analysis in the next section. 

Along the boundary dQ the source potential function 
ipsix) produces the field Hg in the gap interval [—a, -|-a] 
thus Hg = — using the definition eq. (|^. Outside 
[—a, -ba] interval, the potential is constant and the field 
is zero. 


In fact the field is created by a current injected in the 
coil (see fig. [^, as given by Ampere law: NI = Hg with 
I the current and N the number of coil turns. 

The resulting source potential over dQ is given by: 


Psix) 


-K) 

X < —a, 


HgX 

—a < X < a, 

(15) 

Vo 

X > a 



Since psix) is a continuous function we infer that Vq = 


Hga = NIa. 


Over the region Q, the evaluation of the potential 
(p{x,y) requires integration of eq. 12 to get: 


(t>{x,y) = ^ 

TT 


(x + a) tan 


f X + a 

y 


— (x — a) tan 


-1 


X — a 


_ y (37 + a) -b 2/2 


2 (x — a) + 2/2 


(16) 


The field components Hx{x,y),Hy{x,y) over Q ob- tained from definition § are drawn versus x/a for a 
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fixed value of y/a in fig. 



x/a 



x/a 


The solution 4)(fc, y) being a superposition of and 

g-Wv^ the acceptable function for y > 0 is thus 

$(A:,y) = $(fc, The solution (j){x,y) being the 

inverse Fourier Transform of $(fc,y) and since the lat¬ 
ter is the product of two Fourier Transforms <i>(fc,0) and 
g-|fe|y^ then by the convolution theorem, (j){x,y) is the 
convolution of the inverse Fourier Transforms of $(fc, 0) 
and 

Using the boundary condition 4>{x,y = 0) = 4)s{x) 
from eq. and using Fourier Transform definition re¬ 
sults in: 




1 

27r 


r+oo 


(j)s{x)e^^^dx 


Moreover recall that: 


(t>s{x) 

(19) 


j-i[e-|fc|y] 


2?/ 

x'^ + y^ 


f{x,y) 


( 20 ) 


since k and x are Fourier Transform pairs. 

Applying the convolution theorem to the inverse 
Fourier Transforms yields the potential as: 


1 f+°° ^ ^ ^ 

(t>s{x)f[x-X ,y)dx 


' —oo 
^+oo 


1 


- % - dx' (21) 

(a: — x')^ -I- y2 


FIG. 3: Variation of the fields Hx{x,y), Hy{x,y) accord¬ 
ing to the Karlqvist approximation with x for fixed y/a = 
10~^, 10“^, 10“^. The edge effect is not visible in H^{x,y) 
whereas it is weaker than reality for Hy{x,y) (see Conformal 
Mapping section). 


which is exactly the solution eq. 12 found previously 
with 2D Green’s function. 


III. EXACT FOURIER ANALYSIS 


B. Fourier Transform solution 

We retrieve the Green function solution of the 
Boundary-value problem eq.|^by introducing the Fourier 
transform: 

1 /'+°° 

<i>{k,y) = T[<j3{x,y)] = — J (j){x,y)e''’^^dx (17) 

Laplace equation is transformed into: 

^ - k'^^ik,y) = 0 for y > 0 (18) 

oy^ 


Karlqvist solution of BVPl eq. is now reexamined 
such that the gap field is no longer imposed. Thus we 
have another BVP (called BVP2) corresponding to the 
exact solution given by: 

A(/)(r) = 0 in Regions I and II, 

4>{r) = ±Vb on pole surfaces (BVP2) (22) 

The linear gap potential (/)/(a:, y) found previously is 
corrected by adding a Fourier surn^ with unknown dis¬ 
crete coefficients A„ in order to satisfy BVP2 eq. The 
potential (j)ii{x,y) in the upper half-plane is expressed 
with continuous coefficients C{k) obtaining: 
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4>i{x, y) = HgX + ^ An sin(n 7 ra;/a)e"'’^^/“ 0 < x < a, y < O', 

n—0 

poo 

4 >ii{x,y) = / C'(fca) sin(fcx) 0 <x<oo, y >0 

Jo 


(23) 


In order to relate both sets of coefficients An and 
C{ka), we match the potential and derivative at y = 0 
obtaining: 


C{ka) 


I sin(fca) 
Tr {kaf 


n—1 


sin(fca) 

{ka)^ — (n7r)^ 
(24) 



Region I 


FIG. 4: (Color on-line) Regions I and II where A(j>{r) = 0 are 
bounded by red borderlines defining BVP2 eq. with pole 
surfaces considered as two distinct equipotentials (j) = ±Vo. 


Eliminating coefficients C{ka) allows us to derive rela¬ 
tions among A^ coefficients: 





OO 

''^{ — l)^mrAnImn + HglmO 

.n—1 


(25) 


The coefficients Inm are given in the Appendix. The 
above constitute an infinite Algebraic system of linear 
equations that can be solved after numerical truncation 
and performing LU decomposition technique^ (see Ap¬ 
pendix). However it is preferable to use the exact values 
obtained from the combination of Fourier results and con¬ 
formal mapping described in the next section since the 
accuracy of the numerical values of A^ decreases with 
the order m when LU decomposition is performed. 



FIG. 5: Fourier reconstructed potential (j)i{x,y = 0) and field 
Hx{x,y — 0) in the gap. Field oscillations due to Gibbs phe¬ 
nomenon! arise from finite sum of Fourier coefficients. They 
are damped with Sigma Lanczos factors affecting as well the 
edge effect. 


IV. CONFORMAL MAPPING 

The Schwarz-Christoffel transformatioiP is used to 
map the semi-infinite pole region of the z-plane onto the 
real axis in the complex w plane (see fig.[^. 

A polygon with vertices located in the z plane with 
interior angles at is transformed into a sequence of points 
wi,W 2 ,--- along the real axis in the w plane with (see 

fig.§: 

aw 

(26) 

The inspection of polygon ABODE in the z plane (see 
%• [ 2 ]) maps BCD vertices with corresponding angles: 
ai = \,OL 2 = §, 0:3 = ^ into B’C’D’ with abscissae 
Wi = W 2 = 0 , 1 x 3 = 1 with Schwarz-Christoffel relation: 


dz y/w — 1 
dw ^ w 


(27) 


where 7 is a complex constant to be determined. 
Integration of eq. [27| yields: 


z = 2 "f{\/w — 1 — tan ^ \/w — 1 ) + zq (28) 
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FIG. 6 : Schwarz-Christoffel mapping of the polygonal region 
ABCDEF of the jz-plane onto wi,W 2 , ...we on the real axis of 
the complex w plane. 



we use the Schwarz-Christoffel transformation 27 to get: 


[ dz= [ 7 ^^ ^ dw (30) 

Jb Jb' w 

This contour integral is evaluated with the replacement 
w = re*® with r —>■ 0 such that B’ corresponds to 0 = tt 
and C’ to 0 = 0: 


a = lim / 7 


° \/[re*®- 1] 


ire^^dd = ^tt 


(31) 


Thus ') = ajiT and zq = a. 

After conformal mapping, we define an analytic com¬ 
plex functiorPF(z(;) = tp + itj) = lnw-|-zVo built with 
the magnetostatic flux ip and the potential (p satisfying 
the boundary conditions in the w plane, namely: (p = 0 
along A’B’ (since arg(iy) = tt), (p = Vq along C’D’ and 
D’E’ (arg(w) = 0 ) (see fig.[|). 

The magnetic held is obtained from F(w) after differ¬ 
entiating with the complex operators gy = (gy + ^) 
and g| = *(gy — ^) that transform field expression ® 

into a 2D complex form with = —Im 


dF(w) 

dz 


and 


Hy = —Re 


dF{w) 


dz 


The fields Hx{x, y), Hy{x, y) are obtained after relating 
w to z = X + iy using relatio n and drawn versus x/a 
for a fixed value of y/a in fig. S] 

Hx{x, y) and Hy{x, y) behave around the edges accord¬ 
ing to the Sommerfeld radiation conditioiP. We find that 
y“3 for X « ±a which is expected from the 




V 


edge rule stating that the fields near an edge behave as 
with V = —N where p is the distance to the 

edge and a the edge angle. Applied to our case, we have 
y = p,a = f yielding v = | thus the y~ s behaviour 


confirmed by Matzner et al. 


rf 


A. Exact Fourier Am coefficients 


FIG. 7: The polygonal region ABGDE in the 2 plane is 
mapped onto A’B’G’D’E’ along the real axis of the complex 
w plane. Line BC is mapped to the line B’G’ deformed into 
a half-circle with vanishing radius r —>■ 0 to allow reveal the 
angle 0 as 0 for C’ and tt for B’. The potential cp = 0 over 
A’B’ and (p = Vo over G’D’E’. 

Using tan“^ z = —i tanh“^(f 2 :) and tanh“^(iz) = 
iln(^) gives: 


Conformal mapping is now used to extract the ex- 
actP An coefficients that are used in the Fourier solution 
eq.[^ 

Since the gap solution given by coefficients that 
are related to the half-space solution coefficients given 
by C{ka), it might be useful to turn to Fourier transform 
of the fields obtained from the conformal transformation. 

The Fourier transform (cf. definition [l7|) of the head 
magnetic field J^[Hx{x,y)] is related to the potential 
(p{x,y) Fourier transform via C{ka) such thafP^ 


z = 27 


\/w — 1 + 2 


1 -b i^w — 1 


zo 


(29) 


^ 1 — J 

In order to evaluate the unknowns 7 and Zq we use the 
map eq. I^D to D’ thus ic = 1 corresponds to z = a and 


J^[(p{x,y)] = — —C{ka)e'^y, 
T[H,(a:,y)] = ~Cika)e'^y. 


( 32 ) 
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Using ka = rmr in the C{ka) expression 24 yields: 


= —(33) 

mTT 

Taking y = 0 gives 0)] and using 

relation J'[iJa;(a:, 0)] = VoS{ka) with S given by Wilton 
et alM 


S{rmr) = Im 


■j_ r+” s 

2Tr 



(34) 

it suffices to calculate the above integral in order to 
evaluate the exact Am coefficients. 

The integral is straightforwardly calculated by complex 
contour integratiorP. There is a single pole s = i of order 
m + 1 in the upper half-plane. Integrating over the real 
axis and closing the contour by a half-circle in the upper 
half-plane, we obtain from the coefficient of the 

S'(m7r) Laurent series the value with: 


m— 1 


Km - ^ 


n—0 


m — 1\ m + n 


n J (m — n)\ 




(35) 


The coefficients Am = given in the Ap¬ 

pendix and compared to Fourier results. 


FIG. 8: Variation of the fields Hx{x, y), Hy(x, y) with x/a for 
fixed y/a = 10“^, 10“^, 10~^. Sommerfeld radiation singular¬ 
ities at a; = ±a edges are enhanced as y/a decreases. 



FIG. 9: Scaling with y of the edge fields Hx{a,y), Hy(a,y) 
confirming the radiation condition exponent predictecpl as 
1/3. 


V. CONCLUSION AND PERSPECTIVES 

The Karlqvist approximation pioneered study of mag¬ 
netic recording and replay heads. It is a very important 
landmark in the history of recording. 

Karlqvist made a very elaborate approach based on 
a specific BVP (called BVPI) he solved with Green’s 
function techniques that were very popular at the time 
of his work. 

While Karlqvist solution is based on a constant gap 
field, Fourier methods and conformal mapping methods 
tackle the Laplace problem directly as a boundary value 
problem called BVP2 without considering the constant 
gap field assumption to obtain a more accurate and real¬ 
istic description of the head field satisfying Sommerfeld 
edge condition at angular points. 

While Green’s function is important and yields several 
fundamental results, it lacks the flexibility and versatil¬ 
ity of conformal mapping that allows to derive the scal¬ 
ing law at the edge points as well as evaluate the exact 
Fourier expansion coefficients. 

It is remarkable that it took 40 year^ to calculate the 
exact coefficients for a straightforward Laplace problem 
consisting of two semi-infinite equipotential domains sep¬ 
arated by a gap and facing a free half-plane. Using popu¬ 
lar software packages, be they symbolic (like Mathemat- 
ica or Maple) and Finite-Element Method based such as 




























COMSOL-Multiphysics can approach the problem read¬ 
ily with varying degree of success but miss the physical 
points raised in this work. 

In spite of all developments in magnetic recording the¬ 
ory, Karlqvist approximation remains a major milestone 
in recording physics and is always used as a standard 
gauge with respect to any recording theory. 


Appendix A: Evaluation of the exact Fourier 
Coefficients 


The numerical evaluation of the exact is straight¬ 
forward but the evaluation of the infinite algebraic sys¬ 
tem of equations eq. is tricky since it involves non¬ 
trivial intermediate operations such as infinite summa¬ 
tions, evaluation of the Imn,ImO integrals... In addition, 
the presence of minus signs in the system might induce 
numerical conditioning. 

Since the Am coefficients satisfy the relations (taking 

Hg = lm eq.|^: 



2(_i)m+l 

TT 


oo 

^ ^ ( 1) Ajilmn “t” .ImO 

_n—1 


(Al) 


it is possible to writ^for them a linear system AIx — b 
where the matrix M elements are given by: 


— (-lj'"+"n 7 r/ -I- < 4 'm — n, m — 1,2,3... 

IVlmn-i ij mTlmn^<^Q m ^ n, 71=1,2,3... 

(A2) 

the unknowns Xm = Am and the RHS terms are given 
by bm = TO = 1, 2, 3.... 

The Imn terms are integrals given by: 


Imn. — 


xsm X 


'o — (m7r)^][x^ — 


dx 


(A3) 


can be expressed after partial fraction decomposition 
into simpler integrals of the form: 


dm — 


X sin^ x 1 

-dx = - 


10 [x'^ — (tott) ] 


where: 



X 

(tott)^] 


dx - Lm 


(A4) 


L 


m 



X cos 2a; , 

- Trdx 

[a:^ — (tott)^] 


(A5) 


Integrals Lm can be evaluated with contour integra- 
tioiP or obtained from Gradshteyn-Ryzhik table^^^. As 
a result, terms Imn,Imm and Imo are expressed analyti¬ 


cally as: 


Imn. — 


[ln(TO/7i) — Ci{2rmr) Ci{2mr)] 


277^(71^ — TO^) 


ImO — 


Si{2rmr) 
2TO7r 
1 

2vn?TT'^ 


[lE + ln(2TO7r) — Ci{2rmT)\ 


n ^ m, 
n = m, 


77 = 0 


where 'Je = 0.57721566... is Euler-MascheronP^^ con¬ 
stant and Ci{x), Si{x) are the cosine and sine integralJi^^ 
given respectively by: 


^ ^ f°° cos(f) , r, r \ sin(f) 

Ci{x) = — - -dt, Si{x) = / - -dt (A6) 

Jx t Jo t 

In order to solve the system with LU decompositiorP 
we truncate it to an A x iV system with N taken succes¬ 
sively as N = 100,200...1000. The convergence is moni¬ 
tored by comparing the LU coefficients to the exact coef¬ 
ficients as N varies. After establishment of convergence, 
a value of A = 1000 allows us to assess the accuracy in 
the evaluation of Am as to increases (see Table |T]). 


n 

Exact 

LU decomposition 

Wilton et al.^ 

1 

-8.6157113x10"^ 

-8.6155161x10"® 

-8.6157121x10"® 

2 

2.9150246x10"^ 

2.9149190x10"® 

2.9150244x10"® 

3 

-1.5254218x10"^ 

-1.5253429x10"® 

-1.5254219x10"® 

4 

9.5924977x10"^ 

9.5917108x10"® 

9.5924996x10"® 

5 

-6.6803270x10"® 

-6.6796015x10"® 

-6.6803293x10"® 

6 

4.9651614x10"® 

4.9644802x10"® 

4.9651619x10"® 

7 

-3.8608240x10"® 

-3.8602557x10"® 

-3.8608965x10"® 

8 

3.1035715x10"® 

3.1029664x10"® 

3.1035778x10"® 

9 

-2.5590272x10"® 

-2.5585052x10"® 

-2.5590907x10"® 

10 

2.1529666x10"® 

2.1524425x10"® 

2.1530113x10"® 

11 

-1.8415294x10"® 

-1.8406087x10"® 

-1.8411567x10"® 

12 

1.5957345x10"® 

1.5953275x10"® 

1.5958599x10"® 

13 

-1.3977018x10"® 

-1.3985117x10"® 

-1.3990281x10"® 

14 

1.2364594x10"® 

1.2378943x10"® 

1.2383976x10"® 

15 

-1.1009253x10"® 

-1.1049084x10"® 

-1.1054002x10"® 

16 

9.9083269x10""^ 

9.9341292 xl0""‘ 

9.9389255x10""* 

17 

-8.9483510x10""^ 

-8.9890108x10""^ 

-8.9937093x10"* 

18 

8.2353526x10""^ 

8.1800920x10""^ 

8.1846816x10"* 

19 

-7.4235769 xlO""^ 

-7.4817217x10""^ 

-7.4862316x10"* 

20 

7.0718257x10""^ 

6.8741635 xl0"'‘ 

6.8785821x10"* 


TABLE I: Exact, LU decomposition and Wilton et al. evalu¬ 
ated An coefficient^. As n increases, the discrepancy in the 
coefficients originates from decrease in accuracy of the LU 
decomposition. 
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